%% split ETOPO into countries and into connected countries in Europe
% compute ruggedness measure within 0.25 degree cells

clear
close all
clc

set(0,'DefaultFigureWindowStyle','docked','DefaultFigureVisible','on')

% -----------------------
% DEFINE FOLDER LOCATIONS
folders;

%% load country list
country_list_short;

%% We only need to do this if one of the countries is Brazil or Peru
% load TIFF files 
if any(strcmp(country_names(:, 3),'Brazil')) || any(strcmp(country_names(:, 3),'Peru'))
    tiff_files =   dir(fullfile(datafolder_GFC,'*.tif'));
  
    counter = 1;
    for k = 1:length(tiff_files)
        
        [k length(tiff_files)]
        
        [A, R]  = geotiffread( [datafolder_GFC, tiff_files(k).name] ); 
        
        Ygrid = R.LatitudeLimits(2):-R.CellExtentInLatitude:R.LatitudeLimits(1)+R.CellExtentInLatitude;
    
        % X runs West-East
        Xgrid = R.LongitudeLimits(2):-R.CellExtentInLongitude:R.LongitudeLimits(1)+R.CellExtentInLongitude;
    
        % set cell size, in minutes
        cellsize = 0.5;
        NN = round( cellsize/R.CellExtentInLatitude);  % each 0.5 cell will feature NN^2 cells from the map
    
        % north-west corners of each cell
        Ygrid_coarse = round(max(Ygrid)):-cellsize:min(Ygrid);
        Xgrid_coarse = round(min(Xgrid)):cellsize:max(Xgrid);
    
        clc
    
        for j=1:length(Xgrid_coarse)-1
    
            [ j length(Xgrid_coarse) ]
    
            for i=1:length(Ygrid_coarse)
    
                places(counter).Geometry = 'Polygon';
                places(counter).X = [ Xgrid_coarse(j)          Xgrid_coarse(j) Xgrid_coarse(j)+cellsize Xgrid_coarse(j)+cellsize NaN ];
                places(counter).Y = [ Ygrid_coarse(i)-cellsize Ygrid_coarse(i) Ygrid_coarse(i)          Ygrid_coarse(i)-cellsize NaN ];
    
                places(counter).centerX = mean( places(counter).X( ~isnan( places(counter).X ) ) );
                places(counter).centerY = mean( places(counter).Y( ~isnan( places(counter).Y ) ) );
    
                temp = A( NN*( i-1 )+1:NN*i,...
                          NN*( j-1 )+1:NN*j );
    
                places(counter).tree_cover = mean( temp(:) );             % average altitude
                counter = counter+1;
    
            end
        end
    
    end
end

%% retain the section of European grid corresponding to each country
for country_n=1:Ncountries
    
    country_icc = char( countries(country_n) );  % country ICC code
    country = icc2name( country_icc );           % country name

    if strcmp(country, 'Peru') || strcmp(country, 'Brazil') 
        
        % load country_bounds
        if ~ispc()   %% we are in a mac

            % folder to save the outcome
            datafolder_LAC_country =  [ datafolder_LAC,country_icc,'/'];

        else   %% we are in a pc   

            % folder to save the outcome
            datafolder_LAC_country =  [ datafolder_LAC,country_icc,'\' ];

        end

        clc
        disp( ['Country: ',country] )     

        % load bounds 
        file  = [ datafolder_LAC_country,'country_bounds.mat'];
        load(file);

        % keep places within the defined boundaries
        keep_place = inpolygon( cell2mat({places.centerX}'), cell2mat({places.centerY}'),...
                              [ min(country_bounds.X)-1 max(country_bounds.X)+1 ],...
                              [ min(country_bounds.Y)-1 max(country_bounds.Y)+1 ] );
        gfc = places( keep_place );

        % save
        file  = [ datafolder_LAC_country,'gfc.mat'];
        save(file,'gfc');

        % plot

        h=figure;
        mapshow(gfc,'FaceColor','white')
        mapshow(country_bounds,'FaceColor','none')
        title(country)
        if ismac
            saveas(h, [datafolder_GFC,'maps/',country,'_map_gfc.pdf']) 
        else
            saveas(h, [datafolder_GFC,'maps\',country,'_map_gfc.pdf']) 
        end
    
    end
       
end


                    